library(ggplot2)
library(ggrepel)
library(multiplot)

#TWAS data 只取random1做部分结果的展示

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201204")
ASD=read.csv("TWAS.ASD.random.con.1.statis.csv",header=T)[6,]	#这是阿尔兹海默症的数据,缩写应该为AD
BD=read.csv("TWAS.BD.random.con.1.statis.csv",header=T)[6,]
Depr=read.csv("TWAS.Depr.random.con.1.statis.csv",header=T)[6,]
SCZ.BD=read.csv("TWAS.SCZ.BD.random.con.1.statis.csv",header=T)[6,]
ASD$group.source="Alzheimer"
BD$group.source="BD"
Depr$group.source="Depr"
SCZ.BD$group.source="SCZ.BD"
rt=rbind(ASD,BD,Depr,SCZ.BD)
rt$logpvalue=-log10(rt$P.value)
rt=rt[order(rt$P.value),]
rt$num=4:1
																#散点图好看些

p1=ggplot(rt,aes(logpvalue,OR))+geom_point(aes(size=3,color=logpvalue))+
  scale_color_gradient(low = "#3C5488", high = "#E64B35")+
  theme_light(base_size = 15)+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+
  geom_text_repel(data = rt,aes(label = paste(group.source,overlap.num,sep="\n")),size = 5)+
  theme(panel.grid.minor = element_blank())+scale_x_continuous(breaks=c(seq(0,3.5,0.5)),limits=c(0,3.5))+
  xlab("log10(P-value)")+ theme(legend.position = "none")	#这里去掉图例是为了更好地画图，后期会在Ai里复制一个图例过去
  
  
  #807 psyASH associated gene Toppgene GO富集图
  setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201230.GO.enrich")
file=read.csv("807psyASH.associated.gene.toppgene..txt",head=T,sep="\t")
unique(file$Category)
library(ggplot2)

file1=file[file$Category=="GO: Cellular Component",]
file1$num=dim(file1)[1]:1
file1=data.frame(file1$Name,file1$q.value.FDR.B.H,file1$Hit.Count.in.Query.List,file1$num)
names(file1)=c("name","FDR","counts","num")
p2=ggplot(file1,aes(FDR,num))+geom_point(aes(size=counts,color=FDR))+
    scale_color_gradient(low = "#E64B35", high = "#3C5488")+
    scale_y_continuous(breaks = file1$num,labels = file1$name)+
    scale_x_continuous(limits=c(0, 0.04), breaks=seq(0, 0.04, 0.02))+theme_light(base_size = 15)+
    theme(panel.grid.minor = element_blank())+labs(x="FDR",y="",title = "CC")+ theme(legend.position = "none")	#这里去掉图例是为了更好地画图，后期会在Ai里复制一个图例过去

file1=file[file$Category=="GO: Molecular Function",]
file1$num=dim(file1)[1]:1
file1=data.frame(file1$Name,file1$q.value.FDR.B.H,file1$Hit.Count.in.Query.List,file1$num)
names(file1)=c("name","FDR","counts","num")
#file1[3,1]="           phosphotransferase activity"	#phosphotransferase activity, alcohol group as acceptor
file1[7,1]="                     transferase activity, transferring phosphorus-containing groups"	#transferase activity, transferring phosphorus-containing groups
p3=ggplot(file1,aes(FDR,num))+geom_point(aes(size=counts,color=FDR))+
  scale_color_gradient(low = "#E64B35", high = "#3C5488")+
  scale_y_continuous(breaks = file1$num,labels = file1$name)+
  scale_x_continuous(limits=c(0, 0.05), breaks=seq(0, 0.05, 0.02))+theme_light(base_size = 15)+
  theme(panel.grid.minor = element_blank())+labs(x="FDR",y="",title = "MF")+ theme(legend.position = "none")	#这里去掉图例是为了更好地画图，后期会在Ai里复制一个图例过去

file1=file[file$Category=="GO: Biological Process",][c(1:3,6,11,18,21),]
file1$num=dim(file1)[1]:1
file1=data.frame(file1$Name,file1$q.value.Bonferroni,file1$Hit.Count.in.Query.List,file1$num)
names(file1)=c("name","FDR","counts","num")
p4=ggplot(file1,aes(FDR,num))+geom_point(aes(size=counts,color=FDR))+
  scale_color_gradient(low = "#E64B35", high = "#3C5488")+
  scale_y_continuous(breaks = file1$num,labels = file1$name)+
  theme_light(base_size = 15)+
  theme(panel.grid.minor = element_blank())+labs(x="FDR",y="",title = "BP")+ theme(legend.position = "none")	#这里去掉图例是为了更好地画图，后期会在Ai里复制一个图例过去
  
  
library(ggrepel)# gene set
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
result=read.csv("./20201231.GeneSet/GeneSet_enrichment.807psyASH.csv",header = T)
rt=result[result$con.file=="con_genotype.random1.anno",]	#只取random1做部分结果展示
 rt$FDR=p.adjust(rt$p.value,method = "BH")
 rt$logFDR=-log10(rt$FDR)
 p5=ggplot(rt,aes(logFDR,OR))+geom_point(aes(size=3,color=logFDR))+
     scale_color_gradient(low = "#3C5488", high = "#E64B35")+
     theme_light(base_size = 15)+geom_vline(xintercept = -log10(0.1),color="red",linetype="dashed")+
     geom_text_repel(data = subset(rt,FDR<0.1),aes(label = paste(geneset,case_overlap,sep="\n")),size = 5)+
     theme(panel.grid.minor = element_blank())+scale_x_continuous(breaks=c(seq(0,3.5,0.5)),limits=c(0,3.5))+
     scale_y_continuous(breaks=c(seq(0,3,0.5)),limits=c(0,3))+xlab("log10(FDR)")+ theme(legend.position = "none")	#这里去掉图例是为了更好地画图，后期会在Ai里复制一个图例过去

	 
#	GWAS结果
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210114.GWAS.enrich")
sel=list.files(pattern="con.1")[-c(6)]

ADHD=read.csv("ADHD.GWAS.random.con.1.statis.csv",header = T)[6,]
ASD=read.csv("ASD.GWAS.random.con.1.statis.csv",header = T)[6,]
BD=read.csv("bd.GWAS.random.con.1.statis.csv",header = T)[6,]
Depr=read.csv("Depr.GWAS.random.con.1.statis.csv",header = T)[6,]
scz=read.csv("gwas.pgc3.scz.GWAS.random.con.1.statis.csv",header = T)[6,]

ADHD$group="ADHD"
ASD$group="ASD"
BD$group="BD"
Depr$group="Depr"
scz$group="SCZ"

library(ggrepel)
data=rbind(ADHD,ASD,BD,scz)
data=data[order(data$P.value),]
data$group.s="N.S"
data[data$P.value<0.05,]$group.s="S"

data$num=4:1
names(data)=c("Term","counts",names(data)[3:7])
data$log10pvalue=-log10(data$P.value)

p6=ggplot(data,aes(log10pvalue,OR))+geom_point(aes(size=3,color=log10pvalue))+
  scale_color_gradient(low = "#3C5488", high = "#E64B35")+
  theme_light(base_size = 15)+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+
  geom_text_repel(data = data,aes(label = paste(group,counts,sep="\n")),size = 5)+
  theme(panel.grid.minor = element_blank())+
  xlab("-log10(P-value)")+ theme(legend.position = "none")

layout <- matrix(c(1, 2, 3, 7,5,7,4,5,6), nrow = 3)
multiplot(plotlist=list(p1,p6,p2,p5,p3,p4),layout=layout)